Iowa_Regen_RMarkdown

Introduction

# The intention of this RMarkdown is to serve as a repository for all code related to my dissertation research on the impacts of regenerative agriculture adoption in Iowa. I assembled all code written thus far and create a few placeholders for future code I intend to write. My hope is to send this document to collaborators and committee members as needed, and to have an organized place to access my code in the future. 

#The objectives of this project are as follows: 
  # 1. Quantify management impacts on soil C and N stocks in large-scale row crop systems in the U.S. Midwest – with application to the broader scales
  # 2.Beyond C and N stock differences, my goal is to understand the system-level impacts of regenerative practices, as they are currently defined and used by producers in the Midwest. 

# I am achieving this by: 
# - Evaluating C and N stocks on neighboring farms using conventional and regenerative management systems in bulk soil and density-separated fractions 
# - Comparing on-farm results to those generated on five long-term tillage experiments across Iowa   
# - Characterizing microbial community structure between farms under different management 
# - Conducting system-level (biogeochemical and economic) modeling using management data from farmer collaborators 

# In the fall of 2022, I sampled 22 farms throughout Iowa using one of two management systems: 
# 1. Regenerative - annual cover crops (at least 4 years continuous use) and no-till farming (at least 8 years continuous use) 
# 2. Conventional - annual or biannual tillage and bare winter fallow 

# All farms were growing corn and soybeans in rotation, and were planted to soybeans in 2022.

include_graphics("images/conv_regen.PNG")

Different management practies impact soil structure, biodiversity, and soil C and N stocks.

include_graphics("images/study_design.PNG")

Study design - including locations of farms sampled, method of selecting sampling points, soil depths sampled, and data generated for each sample.

# Load C and N concentration data, includes the re-runs done in October 2023
CN_data <- read.csv("data/CN_data_all_08.24.23_final.csv", fileEncoding="UTF-8-BOM")

# find the average of any duplicates 
C_avg <- aggregate(CN_data$per_C, list(CN_data$Sample_ID), mean)
colnames(C_avg) <- c("Sample_ID", "per_C")
N_avg <- aggregate(CN_data$per_N, list(CN_data$Sample_ID), mean)
colnames(N_avg) <- c("Sample_ID", "per_N")

CN_avg <- left_join(C_avg, N_avg, by = "Sample_ID")

#load sample info and Google Sheets info (contains BD values)
GSheets <- read.csv("data/Iowa_Regen_Data_GSheets.csv", fileEncoding = "UTF-8-BOM")

SampleInfo <- read.csv("data/Iowa_Regen_Sample_Info.csv", fileEncoding = "UTF-8-BOM")

# join the CN data sheet with the BD and sample info sheet

# I'm joining them in this order to remove any missing samples (There were several 60-100 samples that we didnt have)
sample_info <- left_join(GSheets, SampleInfo, by = "Sample_ID")

data1 <- left_join(CN_avg, sample_info, by = "Sample_ID")


# replace all NAs with 0's

#add in carbonate data, where applicable
inorgC <- read.csv("data/inorgC.csv", fileEncoding="UTF-8-BOM")

data2 <- left_join(data1, inorgC, by = "Sample_ID") %>%
  mutate(per_inorgC = coalesce(per_inorgC, 0))

data2$per_OC <- data2$per_C - data2$per_inorgC

#add in root data 
roots <- read.csv("data/Root_LOI.csv")

data3 <- left_join(data2, roots, by = "Sample_ID")

#rename columns with weird long names
data3 <- data3 %>% rename(BD_fine = BD.of.fine.soil.fraction..with.rock.correction.)

#add in texture data 
texture <- read.csv("data/Iowa_Regen_Texture.csv")

#rename columns, change #DIV/0's to NAs (indicates analyses are not complete yet)
texture <- texture %>% mutate(Sample_ID = Sample.ID, 
                              Sand_pct = Sand.., 
                              Silt_pct = Silt.., 
                              Clay_pct = Clay..) %>%
 mutate(Clay_pct = ifelse(Clay_pct == "#DIV/0!", NA, Clay_pct))


data4 <- left_join(data3, texture, by = "Sample_ID")

#add in litter data (placeholder, these analyses are not complete yet)

#filter out R03 no cover crops and R14
data5 <- data4 %>% filter(Trt != "badRegen") %>% filter(Trt != "ncRegen")

#add in sheet with simpler pair names 
names <- read.csv("data/Iowa_Regen_new_pair_names.csv")

data6 <- left_join(data5, names, by = "Field_ID")

# Now add in the fractionation data. This is in a slightly different format - it is currently in a long format and I want to transform it to wide.I'm also still generating a lot of data in this category, so I will need to continue reloading and reformatting as needed. 

frac <- read.csv("data/Frac_CN_data.csv", fileEncoding="UTF-8-BOM")

to_remove <- c("_LPOM", "_MAOM", "_HPOM")

frac_select <- frac %>% dplyr::select(Sample_ID, Frac_type, per_N, per_C) %>% mutate(Sample_ID = str_replace_all(Sample_ID, c("_LPOM" = "", "_MAOM" = "", "_HPOM" = "")))

#selecting just the columns I know I will be using

data <- data6 %>% dplyr::select(Sample_ID, Depth_increment, BD_whole_soil, BD_fine, Soil_moist_per, Total_sample_oven_dry_mass, TOTAL_ROCK_mass, Worm_count, Root_picking_method, pH, Field_ID, Pair_num, Region, Trt, Rep_num, Depth, per_C, per_OC, per_N, Root_mass, RootC_g, Sand_pct, Silt_pct, Clay_pct, New_pair_num, New_Field_ID, MAP, MAT)

write.csv(data, "data/Iowa_regen_data_select.csv")

# This part of the script takes a long time, but should only need to be run once because I'm writing rasters for each county. 
# IA_counties <- counties(state = "IA")
# MN_counties <- counties(state = "MN")

# Filter out each county where my sampling points are located. I'm sure there is a way to make this cleaner, but I'm in a hurry
# Adams_co <- IA_counties %>% filter(NAME == "Adams")
# Murray_co <- MN_counties %>% filter(NAME == "Murray")
# Jefferson_co <- IA_counties %>% filter(NAME == "Jefferson")
# Fayette_co <- IA_counties %>% filter(NAME == "Fayette")
# Washington_co <- IA_counties %>% filter(NAME == "Washington")
# Story_co <- IA_counties %>% filter(NAME == "Story")
# Potta_co <- IA_counties %>% filter(NAME == "Pottawattamie")
# Mitchell_co <- IA_counties %>% filter(NAME == "Mitchell")
# Benton_co <- IA_counties %>% filter(NAME == "Benton")
# Jones_co <- IA_counties %>% filter(NAME == "Jones")
# Bremer_co <- IA_counties %>% filter(NAME == "Bremer")

# Pull elevation rasters for each county
# Adams_DEM <- get_elev_raster(Adams_co, z = 14)
# Murray_DEM <- get_elev_raster(Murray_co, z = 14)
# Jefferson_DEM <- get_elev_raster(Jefferson_co, z = 14)
# Fayette_DEM <- get_elev_raster(Fayette_co, z = 14)
# Washington_DEM <- get_elev_raster(Washington_co, z = 14)
# Story_DEM <- get_elev_raster(Story_co, z = 14)
# Potta_DEM <- get_elev_raster(Potta_co, z = 14)
# Mitchell_DEM <- get_elev_raster(Mitchell_co, z = 14)
# Benton_DEM <- get_elev_raster(Benton_co, z = 14)
# Jones_DEM <- get_elev_raster(Jones_co, z = 14)
# Bremer_DEM <- get_elev_raster(Bremer_co, z = 14)

#Write a raster for each county 
# writeRaster(Murray_DEM, "spatial/Murray_co_DEM_14.tif", overwrite = TRUE)
# writeRaster(Jefferson_DEM, "spatial/Jefferson_co_DEM_14.tif", overwrite = TRUE)
# writeRaster(Fayette_DEM, "spatial/Fayette_co_DEM_14.tif", overwrite = TRUE)
# writeRaster(Washington_DEM, "spatial_Washington_co_DEM_14.tif", overwrite = TRUE)
# writeRaster(Story_DEM, "spatial/Story_co_DEM_14.tif", overwrite = TRUE)
# writeRaster(Potta_DEM, "spatial/Potta_co_DEM_14.tif", overwrite = TRUE)
# writeRaster(Mitchell_DEM, "spatial/Mitchell_co_DEM_14.tif", overwrite = TRUE)
# writeRaster(Benton_DEM, "spatial/Benton_co_DEM_14.tif", overwrite = TRUE)
# writeRaster(Jones_DEM, "spatial/Jones_co_DEM_14.tif", overwrite = TRUE)
# writeRaster(Bremer_DEM, "spatial/Bremer_co_DEM_14.tif", overwrite = TRUE)

This next section goes county by county. Since the fields are all different sizes and shapes, I can’t necessarily use the exact same code to crop each county around the points. I need to look at and check each.

#Read in the Adams county raster 
Adams_DEM <- rast("spatial/Adams_co_DEM_14.tif")

#Load in the pair 01 points, which are located in Adams county 
R01_points <- st_read("spatial/r01_points_geojson.geojson")

## Reading layer `r01_points_geojson' from data source 
##   `C:\Users\eellis\OneDrive - Colostate\Enviro_Data_Sci\Iowa_Regen_Bookdown\spatial\r01_points_geojson.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 9 features and 21 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -94.61348 ymin: 40.9365 xmax: -94.60576 ymax: 40.94236
## Geodetic CRS:  WGS 84

C01_points <- st_read("spatial/c01_points_geojson.geojson")

## Reading layer `c01_points_geojson' from data source 
##   `C:\Users\eellis\OneDrive - Colostate\Enviro_Data_Sci\Iowa_Regen_Bookdown\spatial\c01_points_geojson.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 9 features and 21 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -94.60818 ymin: 40.92229 xmax: -94.60538 ymax: 40.92818
## Geodetic CRS:  WGS 84

R01_points_prj <- st_transform(R01_points, st_crs(Adams_DEM))
C01_points_prj <- st_transform(C01_points, st_crs(Adams_DEM))

Site01_points <- rbind(R01_points_prj, C01_points_prj)

st_crs(Site01_points) == st_crs(Adams_DEM)

## [1] TRUE

# Check that the points fall in the raster and that they make sense geographically (do you recognize the topography?)
qtm(Adams_DEM) +
 qtm(Site01_points)

## stars object downsampled to 1302 by 769 cells. See tm_shape manual (argument raster.downsample)

# Create a 8000 (meter?) buffer around the points
Site01_buffer <- st_buffer(Site01_points, dist = 8000)

# Crop the DEM to that buffer to reduce processing time 
Site01_crop <- crop(Adams_DEM, Site01_buffer, mask = TRUE)

Calculating the Topographic Position Index (TPI) at Pair01

# Calculate mean elevation of a neighborhood (window size = 3x3 cells)
mean_elevation <- focal(Site01_crop, w = matrix(1, 3, 3), fun = mean, na.rm = TRUE)

# Calculate Topographic Position Index (TPI)
tpi <- Site01_crop - mean_elevation

# Extract the TPI value at the Site01_points 
tpi_values <- extract(tpi, Site01_points)

# Create a join key to merge Site01_points and the TPI values 
Site01_points <- Site01_points %>% mutate(join = row_number())
tpi_values <- tpi_values %>% mutate(join = row_number()) %>% rename(TPI = Adams_co_DEM_14)

#join the TPI dataframe to the point metadata
TPI_Site01 <- left_join(tpi_values, Site01_points, by = "join")

# Find the mean TPI for each farm 
TPI_R01 <- tpi_values %>%
  filter(ID %in% 1:9) %>%  # Filter for IDs 1 through 9, which correspond to the R01 points
  summarise(mean_value = mean(TPI))

TPI_C01 <- tpi_values %>%
  filter(ID %in% 10:18) %>%  # Filter for IDs 10 through 18, which correspond to the C01 points 
  summarise(mean_value = mean(TPI))

Calculating the average slope of the field

# Crop the DEM to the points for each field 
R01_buffer <- st_buffer(R01_points_prj, dist = 300)
R01_crop <- crop(Adams_DEM, R01_buffer, mask = TRUE)

C01_buffer <- st_buffer(C01_points_prj, dist = 300)
C01_crop <- crop(Adams_DEM, C01_buffer, mask = TRUE)

qtm(C01_crop)

## Use this code to look at the buffers you just created 
# qtm(R01_buffer)+
#   qtm(R01_points_prj)+
#   qtm(C01_buffer) +
#   qtm(C01_points_prj)+
#   qtm(Adams_DEM)


# Obtain slope for the cropped raster 
C01_slope = terrain(C01_crop, v = "slope")
R01_slope = terrain(R01_crop, v = "slope")

# Calculate zonal statistics for the slope raster 
C01_stats = zonal(x = C01_slope, z = C01_crop, fun = "mean")
R01_stats = zonal(x = R01_slope, z = R01_crop, fun = "mean")

#Find the mean slope for the whole field 
C01_mean_slope <- C01_stats %>% summarise(mean = mean(slope, na.rm = TRUE))
R01_mean_slope <- R01_stats %>% summarise(mean = mean(slope, na.rm = TRUE))

#Read in the Benton county raster 
Benton_DEM <- rast("spatial/Benton_co_DEM_14.tif")

#Load in the pair 02 points, which are located in Benton county 
R02_points <- st_read("spatial/r02_points_geojson.geojson")

## Reading layer `r02_points_geojson' from data source 
##   `C:\Users\eellis\OneDrive - Colostate\Enviro_Data_Sci\Iowa_Regen_Bookdown\spatial\r02_points_geojson.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 9 features and 21 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -91.81799 ymin: 41.97535 xmax: -91.81518 ymax: 41.97771
## Geodetic CRS:  WGS 84

C02_points <- st_read("spatial/c02_points_geojson.geojson")

## Reading layer `c02_points_geojson' from data source 
##   `C:\Users\eellis\OneDrive - Colostate\Enviro_Data_Sci\Iowa_Regen_Bookdown\spatial\c02_points_geojson.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 9 features and 21 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -91.81324 ymin: 41.97618 xmax: -91.80878 ymax: 41.97792
## Geodetic CRS:  WGS 84

R02_points_prj <- st_transform(R02_points, st_crs(Benton_DEM))
C02_points_prj <- st_transform(C02_points, st_crs(Benton_DEM))

Site02_points <- rbind(R02_points_prj, C02_points_prj)

st_crs(Site02_points) == st_crs(Benton_DEM)

## [1] TRUE

# Create a 8000 (meter?) buffer around the points
Site02_buffer <- st_buffer(Site02_points, dist = 4000)

# Crop the DEM to that buffer to reduce processing time 
Site02_crop <- crop(Benton_DEM, Site02_buffer, mask = TRUE)

# Check that the points fall in the raster and that they make sense geographically (do you recognize the topography?)
qtm(Site02_crop) +
 qtm(Site02_points)

## stars object downsampled to 926 by 1080 cells. See tm_shape manual (argument raster.downsample)

Calculating the Topographic Position Index (TPI) at Pair02

# Calculate mean elevation of a neighborhood (window size = 3x3 cells)
mean_elevation <- focal(Site02_crop, w = matrix(1, 3, 3), fun = mean, na.rm = TRUE)

# Calculate Topographic Position Index (TPI)
tpi <- Site02_crop - mean_elevation

# Extract the TPI value at the Site02_points 
tpi_values <- extract(tpi, Site02_points)

# Create a join key to merge Site02_points and the TPI values 
Site02_points <- Site02_points %>% mutate(join = row_number())
tpi_values <- tpi_values %>% mutate(join = row_number()) %>% rename(TPI = Benton_co_DEM_14)

#join the TPI dataframe to the point metadata
TPI_Site02 <- left_join(tpi_values, Site02_points, by = "join")

# Find the mean TPI for each farm 
TPI_R02 <- tpi_values %>%
  filter(ID %in% 1:9) %>%  # Filter for IDs 1 through 9, which correspond to the R02 points
  summarise(mean_value = mean(TPI))

TPI_C02 <- tpi_values %>%
  filter(ID %in% 10:18) %>%  # Filter for IDs 10 through 18, which correspond to the C02 points 
  summarise(mean_value = mean(TPI))

Calculating the average slope of the field

# Crop the DEM to the points for each field 
R02_buffer <- st_buffer(R02_points_prj, dist = 90)
R02_crop <- crop(Benton_DEM, R02_buffer, mask = TRUE)

C02_buffer <- st_buffer(C02_points_prj, dist = 90)
C02_crop <- crop(Benton_DEM, C02_buffer, mask = TRUE)

# Check the map to make sure they aren't overlapping and that they account for just the field of interest. Adjust the "dist" as needed.  
qtm(R02_crop)+
qtm(C02_crop)

## Use this code to look at the buffers you just created 
# qtm(R02_buffer)+
#   qtm(R02_points_prj)+
#   qtm(C02_buffer) +
#   qtm(C02_points_prj)+
#   qtm(Benton_DEM)


# Obtain slope for the cropped raster 
C02_slope = terrain(C02_crop, v = "slope")
R02_slope = terrain(R02_crop, v = "slope")

# Calculate zonal statistics for the slope raster 
C02_stats = zonal(x = C02_slope, z = C02_crop, fun = "mean")
R02_stats = zonal(x = R02_slope, z = R02_crop, fun = "mean")

#Find the mean slope for the whole field 
C02_mean_slope <- C02_stats %>% summarise(mean = mean(slope, na.rm = TRUE))
R02_mean_slope <- R02_stats %>% summarise(mean = mean(slope, na.rm = TRUE))

#Read in the Benton county raster 
Fayette_DEM <- rast("spatial/Fayette_co_DEM_14.tif")

#Load in the pair 10 points, which are located in Fayette county 
R10_points <- st_read("spatial/r10_points_geojson.geojson")

## Reading layer `r10_points_geojson' from data source 
##   `C:\Users\eellis\OneDrive - Colostate\Enviro_Data_Sci\Iowa_Regen_Bookdown\spatial\r10_points_geojson.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 9 features and 21 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -91.83367 ymin: 42.94761 xmax: -91.82646 ymax: 42.94925
## Geodetic CRS:  WGS 84

C10_points <- st_read("spatial/c10_points_geojson.geojson")

## Reading layer `c10_points_geojson' from data source 
##   `C:\Users\eellis\OneDrive - Colostate\Enviro_Data_Sci\Iowa_Regen_Bookdown\spatial\c10_points_geojson.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 9 features and 21 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -91.82886 ymin: 42.95018 xmax: -91.8237 ymax: 42.95403
## Geodetic CRS:  WGS 84

R10_points_prj <- st_transform(R10_points, st_crs(Fayette_DEM))
C10_points_prj <- st_transform(C10_points, st_crs(Fayette_DEM))

Site10_points <- rbind(R10_points_prj, C10_points_prj)

st_crs(Site10_points) == st_crs(Fayette_DEM)

## [1] TRUE

# Create a 8000 (meter?) buffer around the points
Site10_buffer <- st_buffer(Site10_points, dist = 2000)

# Crop the DEM to that buffer to reduce processing time 
Site10_crop <- crop(Fayette_DEM, Site10_buffer, mask = TRUE)

# Check that the points fall in the raster and that they make sense geographically (do you recognize the topography?)
qtm(Site10_crop) +
 qtm(Site10_points)

## stars object downsampled to 1179 by 848 cells. See tm_shape manual (argument raster.downsample)

# Calculate mean elevation of a neighborhood (window size = 3x3 cells)
mean_elevation <- focal(Site10_crop, w = matrix(1, 3, 3), fun = mean, na.rm = TRUE)

# Calculate Topographic Position Index (TPI)
tpi <- Site10_crop - mean_elevation

# Extract the TPI value at the Site10_points 
tpi_values <- extract(tpi, Site10_points)

# Create a join key to merge Site10_points and the TPI values 
Site10_points <- Site10_points %>% mutate(join = row_number())
tpi_values <- tpi_values %>% mutate(join = row_number()) %>% rename(TPI = Fayette_co_DEM_14)

#join the TPI dataframe to the point metadata
TPI_Site10 <- left_join(tpi_values, Site10_points, by = "join")

# Find the mean TPI for each farm 
TPI_R10 <- tpi_values %>%
  filter(ID %in% 1:9) %>%  # Filter for IDs 1 through 9, which correspond to the R10 points
  summarise(mean_value = mean(TPI))

TPI_C10 <- tpi_values %>%
  filter(ID %in% 10:18) %>%  # Filter for IDs 10 through 18, which correspond to the C10 points 
  summarise(mean_value = mean(TPI))

# Crop the DEM to the points for each field 
R10_buffer <- st_buffer(R10_points_prj, dist = 100)
R10_crop <- crop(Fayette_DEM, R10_buffer, mask = TRUE)

C10_buffer <- st_buffer(C10_points_prj, dist = 130)
C10_crop <- crop(Fayette_DEM, C10_buffer, mask = TRUE)

# Check the map to make sure they aren't overlapping and that they account for just the field of interest. Adjust the "dist" as needed.  
qtm(R10_crop)+
qtm(C10_crop)

# Obtain slope for the cropped raster 
C10_slope = terrain(C10_crop, v = "slope")
R10_slope = terrain(R10_crop, v = "slope")

# Calculate zonal statistics for the slope raster 
C10_stats = zonal(x = C10_slope, z = C10_crop, fun = "mean")
R10_stats = zonal(x = R10_slope, z = R10_crop, fun = "mean")

#Find the mean slope for the whole field 
C10_mean_slope <- C10_stats %>% summarise(mean = mean(slope, na.rm = TRUE))
R10_mean_slope <- R10_stats %>% summarise(mean = mean(slope, na.rm = TRUE))